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Abstract 

The  fundamental  problem  of  collecting  data  in  the  “best  way”  in  order  to  assure  statistically  efficient 
estimation  of  parameters  is  known  as  Optimal  Experimental  Design.  Many  inverse  problems  consist  in 
selecting  best  parameter  values  of  a  given  mathematical  model  based  on  fits  to  measured  data.  These  are 
usually  formulated  as  optimization  problems  and  the  accuracy  of  their  solutions  depends  not  only  on  the 
chosen  optimization  scheme  but  also  on  the  given  data.  We  consider  an  electromagnetic  interrogation 
problem,  specifically  one  arising  in  an  electroencephalography  (EEG)  problem,  of  finding  optimal  number 
and  locations  for  sensors  for  source  identification  in  a  3D  unit  sphere  from  data  on  its  boundary.  In  this 
effort  we  compare  the  use  of  the  classical  D-optimal  criterion  for  observation  points  as  opposed  to  that 
for  a  uniform  observation  mesh.  We  consider  location  and  best  number  of  sensors  and  report  results 
based  on  statistical  uncertainty  analysis  of  the  resulting  estimated  parameters. 


Keywords:  Electromagnetic  inverse  problems,  optimal  design  in  3D  EEG  analysis,  parameter  estima¬ 
tion,  asymptotic  error  analysis. 
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1  Introduction 


In  a  series  of  recent  works  [6,  7,  8,  11,  12]  several  authors  have  developed  a  design  framework  based  on  the 
Fisher  Information  Matrix  (FIM)  for  a  system  of  differential  equations  to  determine  when  and  where  an 
experimenter  should  take  samples  and  what  variables  to  measure  in  collecting  information  on  a  physical  or 
biological  process  that  is  modeled  by  a  vector  dynamical  system.  This  framework  has  also  been  proposed 
for  use  in  inverse  problem  methodologies  in  the  context  of  dynamical  system  or  mathematical  model  param¬ 
eter  estimation  when  one  investigates  the  sufficiency  of  the  number  of  observations  of  one  or  more  states 
(variables).  Experimental  design  using  the  Fisher  Information  Matrix  (FIM),  which  is  based  on  sensitivity 
functions  (traditional  and  generalized),  is  described  in  [7]  for  the  case  of  scalar  data.  In  [8],  the  authors 
develop  an  experimental  design  theory  using  the  FIM  to  identify  optimal  sampling  times  for  experiments  on 
physical  processes  (modeled  by  an  ODE  system)  in  which  scalar  or  vector  data  will  be  taken.  The  methodol¬ 
ogy  can  be  readily  applied  to  problems  involving  ordinary,  partial  and  delay  differential  equations  dynamics 
but  requires  both  a  mathematical  model  and  a  statistical  model.  Early  efforts  published  in  this  area  were 
concerned  with  parameter  estimation  for  one  dimensional  dynamic  systems  in  one  space  or  time  dimension. 
More  recently,  these  ideas  were  successfully  applied  in  [11]  for  an  experimentally  validated  six-compartment 
HIV  model  and  a  thirty-eight  dimensional  enzyme  kinetics  model  of  the  Calvin  Cycle  in  spinach. 

In  [13]  and  [14]  proof-of-concept  numerical  results  for  a  distributed  parameter  system  in  a  3D  one  layer 
spherical  domain  are  presented  for  several  different  design  criteria  (D-optimal,  SE-optimal,  IGSF-optimal). 
In  this  present  effort,  a  more  general  case  in  3D  is  studied.  Motivated  again  by  classical  problems  in  EEG 
analysis,  we  consider  a  stationary  process  modeled  by  a  Poisson  type  equation  with  multiple  interfaces  given 

by  ~ 


V  •  (k(x)\7u(x,  9))  =  g(x,9)  (1) 

Here  9  £  is  the  parameter  to  be  estimated,  fl  is  the  multilayered  sphere  in  M3,  n(x)  is  a  piecewise  positive 
constant  conductivity  function,  u(x,  9)  G  R,  x  £  dfl  is  the  output  and  g  :  R3+p  — »  R  is  the  source.  For  this 
investigation  we  suppose  that  there  exists  a  real  value  90  such  that  the  equation  (1)  accurately  describes  the 
process. 

For  the  estimation  process,  we  formulate  a  statistical  model  [35]  of  the  form  (absolute  error) 

Y(x)  =  u(x,9o)  +  £(x),  x  €  dfl,  (2) 

where  f?o  is  the  vector  of  true  values  of  the  unknown  parameters  and  £  is  a  vector  random  process  that 
represents  observation  error  for  the  measured  variables.  Realizations  of  the  statistical  model  (2)  can  be 
written  as 


y( x)  =  u(x,  9q)  +  e(x),  x  £  dtd.  (3) 

In  almost  every  real  problem  only  a  discrete  set  of  output  data  is  available.  We  denote  {y\,  ...,yn}  the 
measurements  at  A  =  {x\, . . .  ,xn}  C  dCt.  In  this  context,  the  parameter  value  9q  may  be  estimated  by  an 
Ordinary  Least  Square  (OLS)  procedure  yielding  an  estimate  9.  That  is, 

0  =  argrnin  J\{9) 

8^A 

where  A  is  the  set  of  admisible  parameter  values  and  Ja{9)  denotes  the  square  errors  between  the  measured 
and  simulated  outputs  at  the  observation  points,  i.e. , 


M°)  =  '52\u(xi,8)  ~  Vj\2- 

J=i 


(4) 
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Different  choices  of  the  points  x\,. . .  ,xn  £  A  may  lead  to  different  estimates,  and  thus  it  is  important  to 
look  for  the  optimal  set  of  observation  points  that  will  lead  to  an  accurate  parameter  estimation.  This  is  the 
purpose  of  the  so-called  Optimal  Design  methods,  like  D-Optimal,  SE-Optimal,  E-optimal  design  methods 
(see  [6,  7]).  For  the  sake  of  clarity,  in  the  remainder  of  this  article  we  omit  the  subindexes  A  on  J . 

The  results  presented  here  will  clearly  illustrate  the  advantages  of  performing  some  type  of  design  such 
as  D-Optimal  design  in  connection  to  the  estimation  of  the  parameters  in  the  problem  described  by  (1). 

2  The  EEG  Inverse  Problem 

The  electrical  activity  of  the  brain  consists  of  currents  generated  by  biochemical  sources  at  the  cellular  level. 
The  electric  and  magnetic  fields  that  they  produce  can  be  estimated  by  means  of  Maxwell  equations  (see 
[27,  33]).  Based  on  the  properties  of  the  tissues  involved,  the  velocity  of  propagation  of  the  electromagnetic 
waves  caused  by  potential  changes  within  the  brain  is  such  that  the  effect  may  be  detected  simultaneously  at 
any  point  in  the  brain  or  in  the  surrounding  tissues.  This  fact  justifies  the  use  of  a  static  approximation  of 
Maxwell  equations.  This  approximation  uncouples  the  equations  for  the  magnetic  and  electric  fields  leading 
to  a  3D  Poisson-type  equation  with  interfaces  that  relates  the  electric  potential  u  in  the  head  with  the 
impressed  current  C. 

The  equation  that  relates  the  electric  potential  u  and  the  impressed  current  C  reads: 

V  •  (k(x)\7 u(x))  =  V  •  C,  x  £  D, 

^  =  o  xedn, 

av 

where  the  volume  fl  represents  the  head,  v  is  the  external  normal  vector  and  n  is  the  conductivity  function. 
The  impressed  current  is  often  represented  by  an  electric  dipole,  C(x)  =  q8(x  —  rq),  where  <5  is  the  Dirac 
distribution,  rq  is  a  fixed  point  in  the  brain  which  represents  the  dipole  location,  and  q  is  the  dipole  moment 
[27].  Observe  that,  since  a  static  approximation  is  considered,  the  values  of  u( x)  do  not  depend  on  time; 
instead  they  correspond  to  a  precise  instant  of  the  underlaying  process. 

A  simplified  geometric  model  is  usually  considered  where  the  head  is  represented  by  a  3D  three-layered 
sphere.  The  nested  subsets  fl3,  fl2\fl3  and  D1\D2  correspond  to  the  brain,  skull  and  scalp,  respectively.  The 
location  of  the  dipole  rq  £  fl3  (see  Fig.  1)  is  in  the  brain.  The  function  n(x)  represents  the  conductivity  of 


Figure  1:  3D  Three-layered  Domain 


the  tissues  involved  and  is  usually  considered  as  a  positive  piecewise  constant  function 


«3 

x  £  Sl3, 

K2 

x  £  n2\n3 

«1 

x  £  D3\D2 

0 

x  ^  n, 
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Assuming  that  the  potential  and  its  normal  derivative  multiplied  by  the  conductivity  function  are  con¬ 
tinuous  across  the  transition  surfaces,  the  following  interface  conditions  are  imposed 


„  r  ,  sdu. 

=  o  Hx)-\ 


=  0  *  =  2,3. 


where  [■] 


denotes  the  difference  between  the  values  of  the  functions  inside  the  brackets  through  the  interface 


surfaces  Si  =  di \,i  =  1,  ..,3  that  surround  the  different  subsets. 

In  this  framework,  the  forward  problem  of  EEG  consists  in  finding  the  electric  potential  u  in  the  head 
for  a  given  current  source  C .  Conversely  the  inverse  problem  of  EEG  consists  in  finding  the  location  rq 
of  the  electric  dipole  and  its  moment  q ,  for  a  given  scalp  potential  u.  The  remaining  parameters  of  the 
model  «i,  «2>  k3>  as  well  as  the  radii  of  the  spheres,  are  supposed  to  be  known.  We  recall  that  the  values  of 
u(x),  x  £  dtt,  do  not  depend  on  time,  since  they  are  collected  at  a  precise  instant. 


EEG  Recordings 


EEG  Recordings 


Figure  2:  Forward  and  inverse  problems  of  EEG 


Electroencephalographic  source  localization  by  noninvasive  techniques  is  an  area  of  interest  in  clinic 
epileptology,  in  particular  concerning  models  of  dipolar  and  distributed  sources  for  the  investigation  of  focal 
epilepsy.  In  this  context  the  scalp  data  corresponding  to  an  instant  of  a  seizure  should  be  considered. 
An  accurate  solution  to  the  inverse  problem  could  provide  useful  information  to  determine  the  location 
of  epileptogenic  zones  within  the  brain  from  EEG  recordings.  In  the  last  decade  several  authors  have 
been  working  in  this  area.  Source  models  were  analyzed  in  [20,  34,  38]  while  forward  and  inverse  problem 
solutions  were  studied  in  [28,  31,  32],  among  others.  In  [16],  the  authors  developed  a  method  for  EEG  source 
localization  based  on  rational  approximation  techniques  in  the  complex  plane  and  they  presented  results 
using  simulated  data.  We  refer  to  the  references  therein  for  a  review  of  the  principal  results  concerning  the 
inverse  problem  in  EEG. 

Since  in  practice  the  scalp  potential  u  is  measured  only  at  a  finite  set  of  points  Xi,  ■  ■  ■  ,  xni  on  the  scalp  where 
the  electrodes  are  placed,  the  inverse  problem  of  EEG  consists  in  recovering  rq  and  q  from  u(x i),  •  •  •  u(xn). 
Challenged  by  this  problem,  we  analyze  the  corresponding  mathematical  model  considering  that  only  a 
finite  set  of  potential  values  are  available.  Parameter  estimation  techniques  and  optimal  design  schemes  are 
performed  to  solve  the  inverse  problem  accurately. 


3  Mathematical  model 

Motivated  by  the  above  formulation  we  consider  the  following  second  order  elliptic  partial  differential  equa¬ 
tion  (PDE) 

V  •  (k(i)Vm( x,  9))  =  V  •  F{x ,  0),  in  f 1, 

3u  (5) 

tt{x,6)  =  0  in  ec«2, 
dv 
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where  Q  is  the  unit  ball  of  R3,  and  k  is  a  positive  piecewise  constant  function  defined  by 


K3 

X  £  f I3, 

K2 

x  £  n2\fi3, 

«1 

X  £  \^2 , 

0 

x  £  Cl, 

(6) 


where  Q|  =  Q  and  Q2,  ^3  are  balls  centered  at  the  origin  with  radii  r2,r$  such  that  0  <  <  r2  <  1. 

The  mathematical  formulation  is  completely  defined  with  the  following  interface  and  boundary  conditions: 


Si 


=  0 


t(x) 


du 

dzk 


=  2,3. 


(7) 


du(x ) 
dr] 


0,  x  £  d G 


(8) 


where  [•]si  denotes  the  difference  between  the  values  of  the  functions  inside  the  brackets  through  the  surface 
Si  =  dfli  and  77  is  the  external  normal  vector.  This  type  of  system  appears  in  a  number  of  applied  problems 
such  electrostatic  interrogation  systems,  medical  imaging,  geophysical  exploration  and  nondestructive  testing, 
etc,  [2,  3,  4,  5,  10,  15,  17,  18,  24,  25,  26,  29,  30], 

Once  again  based  on  the  formulations  discussed  above,  we  consider  a  function  F  of  the  form  F(x,  9)  = 
q5(x  —  rq),  where  5  denotes  the  dirac  distribution.  Then  the  source  V  •  F(x,9)  could  represent  for  example 
an  electric  dipole  with  moment  q  =  (<71,92,93)  G  R3  and  location  rq  =  (rqi,rq2,rq3)  £  H3.  The  parameter 
of  the  model  is  then  6  =  ( rq,q )  £  R6.  Existence  and  uniqueness  of  a  solution  u{x,6)  to  this  case  have  been 
studied  in  [19]  and  its  dependence  with  respect  to  O  and  k  appears  in  [36,  37],  respectively. 

In  [19]  the  author  gives  an  explicit  formula  in  terms  of  a  series,  for  u  and  its  derivatives  with  respect  to 
the  parameters  for  spherical  domains.  This  formula  allows  one  to  compute  u  at  any  point  of  f l  given  a  dipole 
that  can  be  located  anywhere  in  Q.  Since  we  are  only  interested  in  computing  u(x)  for  x  in  the  boundary 
dCt  of  the  domain  and  with  a  dipole  located  in  the  innermost  layer,  i.e.,  rq  £  Q3,  we  will  only  recall  the 
expression  of  u  in  this  situation.  The  formula  for  u  given  in  [19]  is 


47ru(a:)  =  9- 

MKII 


s“c“",  +  iiiiSo)' 


(9) 


where  cosw  =  u^y  ^  is  the  cosine  of  the  angle  formed  by  the  observation  point  x,  the  dipole  location  rq 
and  the  origin,  and  So,  Si  are  given  by 


ro 


^](2n+  l)i?„(r0,  r)P'n  (cos  w), 
ri>  1 


and  Si  =  y^(2n  +  l)R'n(ro,r)P^(cosur).  (10) 

n>  1 


Here  ro  =  ||r9||,  r  =  ||ar|| ,  Pn  is  the  n-th  Legendre  polynomial  and  P'n  its  derivative  which  can  be  computed 
from  the  recursive  formula 


P0{x)  =  1,  Pi(x)  =  x, 

( n  +  l)Pn+1(x)  =  (2 n  +  1  )xPn(x)  -  nP„_i(x), 


n  >  2. 


Finally,  rq  £  f I3  implies  r0  <  r3,  hence  the  definitions  of  Rn  and  R'n  given  in  [19]  reduce  to 


Rn(r0,r)  =  - 


{£73(ro,  0)}i2  a 

A 


1 0  < 


and 


K(ro  ,r)  =  - 


{^3(^0,  0)}22  , 

K3  A 


'  0  1 


(11) 

(12) 
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where  A  =  {Ui(ri,r2)U2(r2,r3)U3(rs,0)}22,  and  the  matrix  Uj(ra,rb )  is  defined  as 


(2n+l  )Uj(ra,rb)  = 


(s)' 


nr  a 
Tb 


0  -i- 
0  ^ 


n(n+l)«j 

rb 


We  can  thus  rewrite  and  as 


so  that 


n  +  1 , 


*-‘(r°’r>  =  -(2n+l)«,.4 

c  _  1  Pn(C0SW)„n-l 


and 


and 


n  +  1  ^  \ 

n(ra+i)Kj  ^  if  ra,  n  G  (rj+i,  rj),  j  =  1, 2 


if  rb  =  0,  j  =  3. 


C  _  1  nPn(C0SW)„n-l 

0  ' 


(13) 

(14) 

(15) 


The  formula  (9)  with  So  and  S±  given  by  the  previous  formula  (15)  allow  us  to  compute  easily  the  solution 
to  (5)  when  rq  G  H3  and  x  G  dCl. 


4  Inverse  problem  formulation 

We  consider  the  inverse  problem  or  parameter  estimation  for  the  mathematical  model  described  by  (5).  It 
consists  in  estimating  the  unknown  parameter  9  =  (■ rq,q )  G  R6  from  the  data  yi,..,yn,  corresponding  to  n 
observation  points  at  x\,..,xn  on  the  boundary  dCl. 

We  choose  the  Ordinary  Least  Square  (OLS)  method  to  calculate  estimates  9  of  90;  that  is,  we  minimize 

n 

\u(xAe)  -  %|2>  OgA,  (16) 

j'=i 

where  A  =  {(r,  q)\r  G  03}.  We  compute  the  solution  u  of  our  model  by  means  of  (9). 

The  OLS-estimator  9  is  then 

9  =  argmin  J(9). 

8eA 

Theoretical  results  concerning  existence  and  ill-posedness  of  the  inverse  problem  associated  with  (5  )-(8) 
have  been  studied  in  [21,  22,  23]. 

As  already  noted,  a  statistical  model  is  necessary  to  study  and  implement  inverse  problem  techniques 
properly.  Here  we  take  the  absolute  error  model  (2)  with  realizations  (3).  In  this  context  the  inverse  problem 
consists  of  the  estimation  of  the  unknown  parameters  9q  from  the  data 

Vi  :=  u(xi,90)  +  ti  ,i  =  l,..,n, 

where  we  assume  that  the  additive  noises  e\, ..,  en  are  independent  realizations  of  a  centered  normal  random 
variable  with  variance  a2 . 

Recall  that  9  is  a  realization  of  a  random  variable  0.  It  can  be  proved  that  under  suitable  hypothesis 
(see  [9,  35]),  0  has  asymptotically  normal  distribution 

0  ~  N(90,  ( <j2F(xi , ..,  xn,  90)y1), 

where  F(x  1,  ..,xn,9)  G  R6x6  is  the  usual  Fisher  Information  Matrix  defined  by 


(17) 


The  partial  derivatives  J^( Xk,0 )  are  the  traditional  sensitivity  functions  that,  assuming  smoothness  on  u, 
quantify  the  variations  in  u  with  respect  to  changes  in  the  jth  component  of  the  parameter  6.  A  precise 
discussion  of  the  hypothesis  and  the  approximations  involved  in  the  above  theory  is  given  in  [9] . 

For  our  model  representation,  the  sensitivities  of  u  with  respect  to  8  can  be  easily  computed  by  directly 
differentiating  (9).  The  computation  of  the  derivative  Vgw  of  u  with  respect  to  the  moment  q  of  the  dipole 
is  immediate  from  (9).  For  the  derivative  A7rqu  the  author  in  [19]  found  that 


47rVr  u(x)  =  -iq-—  cosw  —  —  +  <S2  —  2S3  cosw  +  SlAcosw)2)  +  q  ■  x(Sz  —  —  —  S4  cosw)  } 

A)  L  A)  v  r0  r0  J  \  r0  J  J 


-Si 


So 

ro 


^  /Si  So 

+9i - 1 

1  r0  r0 


cos  w  j  +  2:  |  <7  •  —  ^3 - -  —  S4  cos  uj'J  +  q  ■  XS4 1 


(19) 


where  So  and  Si  are  given  in  (15),  and  S3,  S4  are  defined  as 
S4  =  4  J2(2n+  1)Rn(ro,r)P.^(cosui)  = 


1>  1 


«  3 


i>l 


s,  =  -  E(2n  +  ixtrcrj^co^)  =  --  E 
r°  K 3 


A 


Here  P”  is  the  2nd  derivative  of  Pn  that  can  be  calculated  from  (11),  and  one  can  use  the  expressions  (14) 
of  Rn  and  R'n  to  make  the  simplifications,  and  obtain 


S2  =  ^(2n  +  l)S"(r0,r)Pn(cosa;). 
ri>  1 


Then  using 


we  find 


K(ro  ,r)  = 


n(n  +  1) 


Rn(ro,r) - R'n(r0,r) 

ro 


n(n  -  1)  2 

(2n+l)K3A  0 


S2 


1  ^  n(n-  l)P„(cosw)  n_2 

-3  h  A 


4.1  Optimal  Design 

It  is  often  useful  to  have  some  criteria  to  determine  where  samples  should  be  taken  for  any  type  of  inter¬ 
rogation  problem,  especially  those  that  can  be  expensive  and/or  invasive.  This  is  the  goal  of  the  optimal 
design  techniques:  to  search  for  the  optimal  set  of  observation  points  in  order  to  carry  on  the  estimation 
procedures.  Different  criteria  generally  give  rise  to  different  selections  of  observation  points.  In  most  cases 
optimal  design  methods  for  parameter  estimation  problems  choose  the  sampling  distribution  by  minimizing 
a  specific  cost  function  related  to  the  error  or  to  the  accuracy  in  parameter  estimates.  Data  collected  in  this 
optimal  way  will  lead  to  parameter  estimates  with  increased  accuracy. 

In  view  of  the  asymptotic  distribution  (17)-(18)  it  is  natural  to  choose  the  points  Xi  that  minimize 
F(xi, ..,  xn,  6<o )  in  some  sense  (see  [6,  7]).  A  well-known  and  widely  used  optimal  design  method  is  the 
D-optimal  criteria  that  consists  in  minimizing  det  F(x  1,  ...a;n,0)-1.  Geometrically,  this  corresponds  to  min¬ 
imizing  the  volume  of  the  confidence  ellipsoid  for  the  covariance  matrix  Cov  =  a2F^x. 

With  that  purpose  we  choose  the  set  Ad  =  {xi,d  ■  ■  ■  ,x„,d}  C  dfl  of  n  observation  points  where  the 
measurements  {y\,  ...,  yn}  are  to  be  obtained  by  minimizing  the  function 

G(x  1,  ...xn,  80)  =  det  F(x  1,  ...xn,  6»0)_1,  xlt  ...x„  €  dfl 
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starting  with  some  initial  set  of  points  and  considering  9q  as  an  initial  guess  value  for  the  parameter.  After 
having  selected  the  observation  points,  we  perform  OLS  with  the  optimal  set  Ap  and  the  initial  guess  9 q. 
In  the  next  section  we  present  a  detailed  description  of  the  way  our  numerical  calculations  were  performed 
for  the  problems  of  interest  to  us  in  this  presentation. 


5  Numerical  Preliminaries 

Since  we  consider  simulated  data,  we  are  able  to  study  the  actual  errors  and  analyze  the  performance  of 
Optimal  Design  for  an  optimal  selection  of  data. 

The  values  used  in  the  numerical  simulation  are  given  in  the  following  two  tables.  Table  1  contains  the 
values  for  the  mathematical  model  that  are  used  in  each  numerical  experiments.  They  are  fixed  values  of 
the  model  and  they  are  assumed  to  be  known. 


Layer 

r  i 

Ki 

Oi\n2 

1 

0.33 

0.92 

0.0042 

H3 

0.84 

0.33 

Table  1:  Known  model  values  (used  for  all  simulations). 


Two  different  electric  dipoles  were  chosen  to  illustrate  a  general  behaviour.  The  first  potential  dipole 
(Example  1)  is  defined  by  F(x,y,z )  =  (3, 4,  0)  5((x,  y,  z)  —  (0.3,  0.4,  0))).  It  is  a  parallel  dipole  since  its 
moment  q  =  (3,4,0)  is  parallel  to  the  vector  position  rq  =  (0.3, 0.4, 0).  The  second  electric  dipole  (Example 
2)  is  given  by  F(x,y,z)  =  (2,  —  1, 1)  S((x,  y,  z)  —  (0.3,  0.4, 0))).  Note  that  we  kept  the  dipole  position  while 
its  moment  is  significantly  different.  These  parameter  values  and  their  respective  initial  relative  errors  are 
shown  in  Table  2. 


Potential  Dipole  (Case  1) 

Potential  Dipole  (Case  2) 

Position  rq 

Moment  q 

Position  rq 

Moment  q 

True  Parameter  values  9q 

(0.3,  0.4,  0) 

(3,  4,  0) 

(0.3,  0.4,  0) 

(2,  -1,1) 

Initial  Guess  values  9g 

(0.2,  0.2,  0.1) 

(2,  2,-1) 

(0.1,  0.2, -0.1) 

(0,  1,  -1) 

Initial  Relative  Error 

0.4899 

0.8466 

0.6000 

1.4142 

Table  2:  Parameter  values  for  inverse  problem  experiments 


The  observation  data  is  simulated  as 

yi  =  Ui,+6i,  i  =  l,...,n.  (20) 

where  the  values  ut  =  u(xi ,  9q),  i  =  1,  •  •  •  ,  n  are  computed  as  a  sum  of  20  terms  of  the  form  given  in  (9),  for 
9o  =  (0.3,  0.4,  0,  3,  4,  0)  in  the  first  example  and  9q  =  (0.3,  0.4,  0,  2,  —1,  1)  in  the  second. 

The  observation  noise  e,;  at  each  observation  point  Xi  is  calculated  by  the  Matlab  function  randn  with 
standard  deviations  a  =  0,  0.05,  0.1,  0.15  and  0.2  to  produce  noise  free  data  as  well  as  data  sets  with 
nontrivial  increasing  noise  levels. 

Note  that  the  observation  data  is  taken  for  the  outer  surface  of  the  domain  which  is  considered  as  a  unit 
sphere  (r3  =  1).  Hence,  the  observation  points  can  be  written  as  (x,  y,  z)  =  (sin(a)  cos (</>),  sin(a)  sin(^),  cos(a)), 
with  (a,  (j>)  €  [0, 7r]  x  [0,  2u\.  We  define  a  uniform  grid  on  spherical  coordinates  with  constant  step  Aa  =  7r/30 
for  a  and  A (f>  =  2  *  7r/30  for  </>.  This  gives  us  31  x  31  grid  points  in  the  rectangle  [0, 7r]  x  [0,  27t]  of  the  form 


Figure  3:  Rectangular  grid  in  (</>,  a )  (left)  and  Spherical  grid  (right). 


(a*,  4>j )  where  cc*  =  (i  —  1)  *  Aa,  i  =  1, 31,  4>j  =  (j  —  1)  *  A^>,  j  =  1, 31.  The  gridpoints  are  numbered 
on  the  rectangular  grid  from  left  to  right  and  from  bottom  to  top  as  shown  in  Fig.  3  (left).  The  spherical 
gridpoints  are  shown  in  Fig.  3,  on  the  right.  Note  that  the  bottom  line  of  gridpoints  correspond  to  a  =  0 
and  the  top  line  to  a  =  <f>.  Each  of  these  lines  correspond  to  only  one  point  on  the  unit  sphere,  they  are 
(0,  0, 1)  and  (0,  0,  —1)  in  cartesian  coordinates,  which  are  the  poles.  On  the  other  hand,  the  extreme  points 
of  the  same  line  are  the  same,  due  to  the  periodicity  of  trigonometrical  functions.  One  needs  to  take  these 
facts  into  account  when  considering  observation  points. 

As  mentioned  before,  the  estimated  value  for  the  parameter  is  based  on  the  initial  set  A  of  observation 
points. 

We  ran  numerical  experiments  considering  different  numbers  of  measurements  from  n  =  7  to  n  =  20.  For 
each  case,  the  points  Xi,  i  =  1,  ..,n  in  A  are  chosen  from  the  set  of  observable  gridpoints  as  follows.  First, 
we  fix  the  first  and  the  last  points  of  A  at  the  gridpoints  numbered  as  100  and  650,  respectively;  that  is 
x\  =  xgiQo,  xn  =  xgeso-  In  this  way  we  are  sure  that  they  fall  far  from  the  poles.  The  rest  of  the  points  are 
uniformly  distributed  on  the  numbered  grid,  which  means 

550 

Xi  =  xgk,  k  =  100  +  (l  —  1) - l  =  l,...,n.  (21) 

n  —  1 


6  Parameter  estimation 

We  obtain  two  estimates  9\  and  9d  of  9q  performing  OLS  with  the  initial  guess  9g  and  two  different  sets  of 
observations  points  x\, ...  ,xn: 

•  the  estimate  9\  is  obtained  performing  OLS  with  the  initial  guess  9g  and  the  set  A  =  {ii, . . .  ,xn}  of 
observation  points  given  by  (21). 

•  we  obtain  a  second  estimate  9d  of  9q  applying  OLS  using  the  observation  points  arising  from  the 
D-optimal  design  criteria  as  follows.  We  first  look  for  a  set  An  =  {x® , . . . ,  X®}  of  n  observation  points 
minimizing  the  function 

det  F(x i,  ...xn,  9g)~1 
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starting  with  A  as  initial  observation  points.  We  then  perform  OLS  with  9g  as  initial  guess  for  6  and 
the  observation  points  in  . 

The  D-optimal  set  of  observations  points  are  calculated  by  means  of  the  Matlab  function  fmincon  with 
the  following  optimset:  MaxIter=50000,MaxFunEvals=700000,TolFun=le-21,TolX=le-6,  with  constraints 
0,  7 r  for  a  and  0,  2-7T  for  <f> . 

The  OLS-estimate  is  computed  by  the  Matlab  function  Isqnonlin  with  the  same  options  as  fmincon 
except  that  in  this  case  we  set  TolFun=le-10. 

We  repeat  these  procedures  N  times  (generating  a  new  set  of  perturbations  ei , ,  en  each  time)  thus  we 
obtain  2N  different  estimates:  N  of  them  denoted  by  9A  =  (f^A,  gA),  j  =  1, . . . ,  N  (  OLS  estimates)  and  the 
remaining  ones  denoted  by  93D  =  (rqD,  q0D),  j  =  1, ...  ,N  (  D-Optimal  -OLS  estimates).  We  then  compute 
the  relatives  errors  for  rq  and  q: 


el(r9)  := 


el  (?)  := 


II 9a  -  90 1|  . 


'  90  I 


l9o| 


e3D{rq)  := 


and  then  average  them: 


eA( 


eAl 


' Lf  yu  ii 

IM  ’ 

eD(9)  :=  - 

II  9o  II  ’  ' 

i  N 

3= 1 

eD{rq) 

1  N 

i=i 

(22) 

1  N 

3= 1 

er»(g)  = 

i  N 

=  ^Ee^(9). 

3  =  1 

(23) 

6.1  Example  1:  (parallel  dipole) 

We  first  consider  a  parallel  dipole  where  the  “true”  values  for  the  parameters  are  given  by 

90  =  (0.3,  0.4,  0,  3,  4,  0) 


and  the  initial  “guess”  values  are 

9g  =  (0.2,  0.2,  0.1,  2,  2,  -1). 

Tables  3  and  4  give  the  resulting  values  for  e\ (?’g),  e_o(rg),  and  for  e\(q),  en{q),  respectively,  for  N  =  500 
with  n  =  7, . . . ,  20  observation  points  and  noise  standard  deviations  er  =  0,  0.05,  0.1,  0.15,  0.2  for  Example 
1  (the  parallel  dipole  example).  These  values  are  displayed  graphically  in  Figures  4  (No  optimal  design)  and 
5  (D-optimal  design). 

We  denote  by  9a  =  (fq,\ ,  9a )  and  9d  =  (rqD,q]j)  the  average  of  the  estimated  parameters  {^a}i<1<x 
and  {9jd}i<j<n,  namely 


0a 


i- 1 


Qd 


i—1 


To  compare  the  statistical  relevance  of  9 a  and  9^>  we  compute  the  confidence  intervals  (Cl)  for  6*o  based 
on  these  two  estimates  in  the  following  way.  We  first  generate  two  sets  of  simulated  observations  {ui:D}i<i<n 
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n\a 

cr  =  0.2 

cr  =  0.15 

cr  =  0.1 

cr  =  0.05 

(7  =  0 

n 

eD(rq) 

eA  (rq) 

eD(rq) 

eA  (rq) 

eD(rq) 

eA(r?) 

eD(rq) 

eA(r-g) 

eD(rq) 

eA  (rq) 

7 

0.2855 

1.1511 

0.2053 

1.0705 

0.1299 

0.8931 

0.0654 

0.6581 

0.65e-15 

0.49e-ll 

8 

0.2079 

0.6825 

0.1549 

0.5100 

0.0995 

0.3342 

0.0498 

0.1671 

0.47e-12 

0.76e-ll 

9 

0.2024 

0.4586 

0.1535 

0.3454 

0.1017 

0.2370 

0.0493 

0.1121 

0.16e-14 

0.11e-14 

10 

0.1990 

1.2967 

0.1362 

1.1389 

0.0931 

0.9728 

0.0453 

0.5965 

0.10e-12 

0.61e-12 

11 

0.1946 

0.3172 

0.1453 

0.2392 

0.0958 

0.1481 

0.0467 

0.0728 

0.27e-15 

0.77e-15 

12 

0.1790 

0.3753 

0.1328 

0.2595 

0.0845 

0.1649 

0.0404 

0.0821 

0.74e-15 

0.54e-13 

13 

0.1837 

0.3752 

0.1297 

0.2846 

0.0896 

0.1863 

0.0432 

0.0951 

0.60e-14 

0.20e-12 

14 

0.1636 

0.3295 

0.1213 

0.2175 

0.0793 

0.1375 

0.0400 

0.0722 

0.25e-14 

0.20e-12 

15 

0.1601 

0.2622 

0.1232 

0.1996 

0.0793 

0.1294 

0.0402 

0.0629 

0.14e-14 

0.56e-14 

16 

0.1521 

0.2581 

0.1181 

0.1865 

0.0744 

0.1252 

0.0384 

0.0603 

0.41e-15 

0.77e-15 

17 

0.1546 

0.2990 

0.1153 

0.2133 

0.0732 

0.1385 

0.0371 

0.0676 

0.34e-15 

0.90e-15 

18 

0.1489 

1.0426 

0.1091 

0.9104 

0.0727 

0.7076 

0.0373 

0.4087 

0.16e-15 

0.38e-13 

19 

0.1481 

0.8855 

0.1063 

0.7488 

0.0720 

0.6297 

0.0355 

0.4026 

0.48e-12 

0.13e-10 

20 

0.1360 

0.2067 

0.1019 

0.1477 

0.0695 

0.1024 

0.0350 

0.0494 

0.12e-15 

0.36e-15 

Table  3:  Example  1.  Mean  relative  errors  e£>(rg)  and  e\ (rq)  for  rq  as  defined  in  (22)  for  different  number  of 
observation  points  n  and  different  value  of  the  noise  standard  deviation  cr. 


n\a 

o-  =  0.2 

cr  =  0.15 

o-  =  0.1 

cr  =  0.05 

(7  =  0 

n 

ez>(r9) 

eA(r\j) 

eD{rq) 

eA(r\j) 

eD{rq) 

eA(r\j) 

eD{rq) 

eA(r-g) 

eD{rq) 

eA  (rq) 

7 

0.1151 

0.7379 

0.0861 

0.6432 

0.0591 

0.4887 

0.0289 

0.3058 

0.17e-15 

0.19e-ll 

8 

0.1033 

0.2099 

0.0769 

0.1509 

0.0512 

0.1013 

0.0264 

0.0497 

0.07e-12 

0.21e-ll 

9 

0.1030 

0.1285 

0.0788 

0.0947 

0.0523 

0.0638 

0.0252 

0.0307 

0.02e-14 

0.06e-14 

10 

0.0986 

0.7908 

0.0709 

0.6081 

0.0485 

0.4683 

0.0239 

0.2570 

0.04e-12 

0.51e-12 

11 

0.0921 

0.0825 

0.0741 

0.0677 

0.0481 

0.0425 

0.0235 

0.0216 

0.28e-15 

0.09e-15 

12 

0.0919 

0.0963 

0.0658 

0.0678 

0.0444 

0.0448 

0.0215 

0.0218 

0.16e-15 

0.15e-13 

13 

0.0862 

0.1717 

0.0630 

0.1197 

0.0435 

0.0837 

0.0219 

0.0412 

0.08e-14 

0.12e-12 

14 

0.0807 

0.0834 

0.0617 

0.0586 

0.0400 

0.0387 

0.0207 

0.0192 

0.01e-14 

0.05e-ll 

15 

0.0779 

0.0783 

0.0585 

0.0565 

0.0396 

0.0376 

0.0205 

0.0183 

0. Ole-14 

0.16e-14 

16 

0.0786 

0.0761 

0.0600 

0.0562 

0.0390 

0.0370 

0.0204 

0.0177 

0.09e-15 

0.10e-15 

17 

0.0765 

0.0900 

0.0572 

0.0647 

0.0375 

0.0422 

0.0192 

0.0213 

0.00e-15 

0.40e-15 

18 

0.0744 

0.2039 

0.0571 

0.1697 

0.0374 

0.1323 

0.0185 

0.0714 

0.25e-15 

0.12e-13 

19 

0.0714 

0.5487 

0.0559 

0.4251 

0.0368 

0.3129 

0.0178 

0.1762 

0.15e-12 

0.02e-10 

20 

0.0706 

0.0741 

0.0519 

0.0538 

0.0345 

0.0357 

0.0183 

0.0183 

0.14e-15 

0.03e-15 

Table  4:  Example  1.  Mean  relative  errors  en{q)  and  e\{q)  for  q  as  defined  in  (23)  for  different  number  of 
observation  points  n  and  different  value  of  the  noise  standard  deviation  a. 


and  {wj,A}i<i<n  defined  in  (20)  using  the  set  of  observation  points  Ad  and  A  respectively.  We  then  compute 
the  estimated  variances  a\  and  g2d  defined  as 

^  ti  2  1  ^  2 

o-i=  ,  ajj  =  (ui>D  -u(xP,0D)^  . 

i—  1  2—1 


11 


The  standard  errors  SE(9A ),  SE{9d)  £  R6  are  then  defined  as 

SEl(9A)  =  a\{F(x  i,...,xn;9A)-1)kk,  SEl(0D)  =  a2D(F(x? , . . . ,  x%]  9D)~l)kk 

where  k  =  1, ...  ,6,  and  F  is  the  Fisher  matrix  defined  in  (18).  The  approximate  Cl  at  the  (1  —  a)%  level 
for  the  k-th  component  do,fc  of  do  corresponding  to  9A  and  9e>  are  then  given  respectively  by 

[0A,k  ~  ti-a/2  SEk{9A)  ,  $A,fc  +  tl-a/2  SEk(9A)  ]  (24) 

[0D,k  -  h-a/2  SEk(9o)  ,  9o,k  +  tx-a/2  SEk{9D)}.  (25) 

The  results  for  rq{  1)  (the  first  component  of  rq)  and  q(  1)  (the  first  component  of  q )  are  shown  in  Figures 
(6)-(9).  Figure  6  depicts  the  results  when  a  =  0.05.  From  top  to  bottom:  The  mean  relative  errors,  length 
of  the  confidence  intervals  (built  using  formula  (24)  (in  blue)  and  (25))  and  the  confidence  intervals.  On  the 

left  are  the  results  for  rq(  1)  and  on  the  right  are  the  results  for  q(  1).  Similar  results  for  a  =  0.1,  0.15,  0.2 

are  depicted  in  Figures  7-9. 


Figure  4:  Mean  relative  errors  eA{rq)  and  eA{q)  (no  optimal  design)  as  functions  of  n  and  a  . 


Figure  5:  Mean  relative  errors  ejj(rq)  and  eu(q)  (using  D-Optimal  design)  as  functions  of  n  and  a. 
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Figure  6:  Results  for  noisy  realizations  with  a  =  0.05  using  D-optimal  design  (red)  and  no  optimal  design 
(blue).  From  top  to  bottom:  Mean  relative  errors  en(rq)  and  e\ (rq)  (right)  and  eo{q)  and  e\ (q)  (left); 
length  of  confidence  intervals  for  rq{  1)  (right)  and  g(l)  (left);  confidence  intervals  at  90  %  confidence  level 
for  rq(  1)  (right)  and  q(  1)  (left). 
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Figure  7:  Results  for  noisy  realizations  with  a  =  0.1  using  D-optimal  design  (red)  and  no  optimal  design 
(blue).  From  top  to  bottom:  Mean  relative  errors  e£>(rg)  and  e\{rq)  (right)  and  eo(q)  and  e\ (q)  (left); 
length  of  confidence  intervals  for  rq(  1)  (right)  and  g(l)  (left);  confidence  intervals  at  90  %  confidence  level 
for  rq(  1)  (right)  and  q(  1)  (left). 
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Figure  8:  Results  for  noisy  realizations  with  a  =  0.15  using  D-optimal  design  (red)  and  no  optimal  design 
(blue).  From  top  to  bottom:  Mean  relative  errors  e£>(rg)  and  e,\(rq)  (right)  and  eo(q)  and  e\ (q)  (left); 
length  of  confidence  intervals  for  rq(  1)  (right)  and  g(l)  (left);  confidence  intervals  at  90  %  confidence  level 
for  rq(  1)  (right)  and  q(  1)  (left). 
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Figure  9:  Results  for  noisy  realizations  with  a  =  0.2  using  D-optimal  design  (red)  and  no  optimal  design 
(blue).  From  top  to  bottom:  Mean  relative  errors  e£>(rg)  and  e\{rq)  (right)  and  eo(q)  and  e\ (q)  (left); 
length  of  confidence  intervals  for  rq(  1)  (right)  and  g(l)  (left);  confidence  intervals  at  90  %  confidence  level 
for  rq(  1)  (right)  and  q(  1)  (left). 
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6.2  Example  2:  (nonparallel  dipole) 

To  illustrate  our  findings  with  a  second  example,  we  consider  Example  2  introduced  above  where  the  moment 
and  location  are  non  parallel.  In  this  case  we  chose  true  values  for  do  =  (0.3,  0.4,  0,  2,  — 1,  1)  and  used  initial 
guess  parameter  values  6g  =  (0.1,  0.2, -0.1,  0,  1,  -1).  Results  for  different  levels  of  noise  are  given  in  Tables 
5  and  6.  Findings  for  mean  relative  errors  and  confidence  intervals  of  level  90  %  are  given  in  Figures  12  -  15. 


n\a 

cr  =  0.2 

cr  =  0.15 

cr  =  0.1 

cr  =  0.05 

a  =  0 

n 

eD{rq) 

e\(rq) 

eD{rq ) 

eA  {rq) 

eD{rq) 

eA  {rq) 

eD{rq) 

eA  (rq) 

eD(rq) 

eA  (rq) 

7 

0.6689 

1.6881 

0.5211 

1.6300 

0.3268 

1.6121 

0.1615 

1.4561 

0.77e-015 

1.8753 

8 

0.6009 

1.2706 

0.4581 

1.2127 

0.2975 

1.0428 

0.1480 

0.6399 

0.51e-015 

0.70e-13 

9 

0.4500 

0.6348 

0.3317 

0.4826 

0.2133 

0.3018 

0.1061 

0.1496 

0.81e-015 

0.28e-15 

10 

0.4272 

1.5717 

0.3194 

1.4599 

0.2005 

1.2102 

0.0951 

0.6284 

0.23e-13 

0.56e-13 

11 

0.3774 

0.5114 

0.2665 

0.1378 

0.1790 

0.2730 

0.0886 

0.1257 

0.79e-15 

0.26e-15 

12 

0.4544 

0.4597 

0.3227 

0.3693 

0.2077 

0.2291 

0.1082 

0.1103 

0.12e-14 

0.85e-15 

13 

0.3573 

0.9601 

0.2855 

0.7718 

0.1748 

0.5196 

0.0898 

0.2374 

0.14e-13 

0.47e-10 

14 

0.3666 

0.4256 

0.2572 

0.3180 

0.1653 

0.2064 

0.0833 

0.1004 

0.17e-14 

0.40e-15 

15 

0.3891 

0.4339 

0.2700 

0.3063 

0.1865 

0.2006 

0.0890 

0.0998 

0.73e-12 

0.47e-12 

16 

0.3235 

0.4291 

0.2384 

0.3217 

0.1589 

0.2043 

0.0766 

0.0983 

0.30e-14 

0.65e-15 

17 

0.3107 

0.4009 

0.2277 

0.3129 

0.1559 

0.1928 

0.0758 

0.0988 

0.17e-13 

0.11e-13 

18 

0.3306 

1.079 

0.2494 

0.8558 

0.1645 

0.6944 

0.0779 

0.4231 

0.49e-15 

0.36e-10 

19 

0.3067 

1.225 

0.2317 

1.0808 

0.1562 

0.7549 

0.0783 

0.3312 

0.33e-15 

0.34e-14 

20 

0.2957 

0.444 

0.2294 

0.3061 

0.1429 

0.2042 

0.0719 

0.0976 

0.14e-15 

0.38e-ll 

Table  5:  Example  2.  Mean  relative  errors  eo(q)  and  e\{q)  for  q  as  defined  in  (22)  for  different  number  of 
observation  points  n  and  different  value  of  the  noise  standard  deviation  a. 


n\a 

cr  =  0.2 

cr  =  0.15 

cr  =  0.1 

cr  =  0.05 

(7  =  0 

n 

e-D{q) 

eA(<?) 

eu(g) 

e\(q) 

ec(<z) 

eA(g) 

eu(g) 

eA  (q) 

eD(q) 

eA(«) 

7 

0.2550 

1.4185 

0.1861 

1.3604 

0.1156 

1.3236 

0.0544 

1.1651 

0.25e-15 

1.2367 

8 

0.2319 

0.6808 

0.1727 

0.6256 

0.1073 

0.4962 

0.0513 

0.2662 

0.27e-15 

0.29e-13 

9 

0.2013 

0.2316 

0.1573 

0.1705 

0.0961 

0.1140 

0.0495 

0.0577 

0.15e-15 

0.04e-15 

10 

0.2007 

1.233 

0.1482 

1.0346 

0.0902 

0.8368 

0.0448 

0.3594 

0.04e-13 

0.42e-13 

11 

0.1844 

0.1843 

0.1311 

0.1378 

0.0902 

0.0900 

0.0409 

0.0437 

0.24e-15 

0.22e-15 

12 

0.1920 

0.1851 

0.1351 

0.1435 

0.0851 

0.0912 

0.0447 

0.0459 

0.04e-14 

0.22e-15 

13 

0.1648 

0.5212 

0.1246 

0.4082 

0.0795 

0.2526 

0.0396 

0.1057 

0.05e-13 

0.16e-10 

14 

0.1516 

0.1628 

0.1107 

0.1211 

0.0728 

0.0803 

0.0360 

0.0377 

0.02e-14 

0.00e-15 

15 

0.1644 

0.1713 

0.1202 

0.1244 

0.0806 

0.0824 

0.0388 

0.0410 

0.27e-12 

0.16e-12 

16 

0.1416 

0.1477 

0.1087 

0.1155 

0.0683 

0.0724 

0.0358 

0.0378 

0.13e-14 

0.29e-15 

17 

0.1459 

0.1484 

0.1088 

0.1121 

0.0715 

0.0718 

0.0349 

0.0359 

0. Ole-13 

0.00e-13 

18 

0.1441 

0.322 

0.1045 

0.2434 

0.0680 

0.1924 

0.0333 

0.1198 

0.10e-15 

0.18e-10 

19 

0.1296 

0.944 

0.0993 

0.7791 

0.0656 

0.5600 

0.0344 

0.2457 

0.16e-15 

0.54e-14 

20 

0.1299 

0.197 

0.0984 

0.1377 

0.0641 

0.0917 

0.0321 

0.0441 

0.20e-15 

0.15  e-11 

Table  6:  Example  2.  Mean  relative  errors  en{q)  and  e\{q)  for  q  as  defined  in  (23)  for  different  number  of 
observation  points  n  and  different  value  of  the  noise  standard  deviation  cr. 
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Figure  10:  Mean  relative  errors  e\ (rq)  and  e\ (9)  (no  optimal  design)  as  functions  of  n  and  er. 


Figure  11:  Mean  relative  errors  e£>(rg)  and  en(q)  (using  D-optimal  design)  as  functions  of  n  and  <r. 
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Figure  12:  Results  for  noisy  realizations  with  a  =  0.05  using  D-optimal  design  (red)  and  no  optimal  design 
(blue).  From  top  to  bottom:  Mean  relative  errors  e£>(rg)  and  e\ (rq)  (right)  and  eo(g)  and  e\ (q)  (left); 
length  of  confidence  intervals  for  rq(  1)  (right)  and  g(l)  (left);  confidence  intervals  at  90  %  confidence  level 
for  rq(  1)  (right)  and  q(  1)  (left). 
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Figure  13:  Results  for  noisy  realizations  with  a  =  0.1  using  D-optimal  design  (red)  and  no  optimal  design 
(blue).  From  top  to  bottom:  Mean  relative  errors  e£>(rg)  and  e\ (rq)  (right)  and  eo(q)  and  e\ (q)  (left); 
length  of  confidence  intervals  for  rq(  1)  (right)  and  g(l)  (left);  confidence  intervals  at  90  %  confidence  level 
for  rq(  1)  (right)  and  q(  1)  (left). 
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Figure  14:  Results  for  noisy  realizations  with  a  =  0.15  using  D-optimal  design  (red)  and  no  optimal  design 
(blue).  From  top  to  bottom:  Mean  relative  errors  e£>(rg)  and  e\ (rq)  (right)  and  eo(g)  and  e\ (q)  (left); 
length  of  confidence  intervals  for  rq(  1)  (right)  and  g(l)  (left);  confidence  intervals  at  90  %  confidence  level 
for  rq(  1)  (right)  and  q(  1)  (left). 
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Figure  15:  Results  for  noisy  realizations  with  a  =  0.2  using  D-optimal  design  (red)  and  no  optimal  design 
(blue).  From  top  to  bottom:  Mean  relative  errors  e£>(rg)  and  e\ (rq)  (right)  and  eo(q)  and  e\ (q)  (left); 
length  of  confidence  intervals  for  rq(  1)  (right)  and  g(l)  (left);  confidence  intervals  at  90  %  confidence  level 
for  rq(  1)  (right)  and  q(  1)  (left). 
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7  Conclusions 


In  this  paper  we  have  investigated  a  typical  interrogation  problem  such  as  those  used  in  EEG.  We  have 
shown  the  value  of  using  some  type  of  optimal  design  criterion  (such  as  those  studied  in  [6,  7,  8,  11,  12])  in 
determining  how  to  best  collect  data.  From  the  numerical  results,  we  would  conclude  that: 

•  D-optimal  design  techniques  provide  a  set  observations  points  leading  to  a  more  accurate  estimate  of 
the  parameters  of  interest.  The  results  of  this  paper  emphatically  demonstrate  the  benefits  of  using 
some  type  of  optimal  design  criterion  (D-optimal  in  this  case)  in  deciding  how  data  should  be  collected 
in  a  specific  application. 

•  For  the  specific  EEG  problem  investigated,  the  length  of  the  confidence  intervals  as  well  as  the  mean 
relative  error  do  not  decrease  significantly  for  more  than  10  or  11  observation  points.  Hence  an  optimal 
array  of  sensors  of  this  size  is  sufficient  in  practice.  Moreover,  there  are  dramatic  differences  between 
estimation  accuracies  with  a  smaller  number  of  sensors  (  «  7  or  8)  and  the  optimal  values  of  10  or  11 
sensors. 
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